//
// Sample delta variances of zero set counts. Use hprim.gp to
// process and calculate H'.
//
// Copyright (C) 2009-2011 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#include "common.cpp"

const char * DEFAULT_OUTPUT = "hprim.csv";

static void usage (const char *argv0)
{
  fprintf (stderr, "%s: %s [-H exponent] [-r seed] [-o %s] [-g generator] [-l growthpow startlayer endlayer] [-f fieldsize] [-i interval] [-s samplecount] [-c step]\n", argv0, argv0, DEFAULT_OUTPUT);
  fprintf (stderr, "    -H exponent - set the Hurst exponent\n");
  fprintf (stderr, "    -r seed - set the random seed\n");
  fprintf (stderr, "    -o %s - output data to file\n", DEFAULT_OUTPUT);
  fprintf (stderr, "    -g generator - use the given random generator, put 'list' to see available\n");
  fprintf (stderr, "    -l growthpow startlayer endlayer - specify the layer growth base and layer range\n");
  fprintf (stderr, "    -f fieldsize - specify the field size\n");
  fprintf (stderr, "    -i interval - field sample interval\n");
  fprintf (stderr, "    -s samplecount - sampling count\n");
  fprintf (stderr, "    -c step - discretion step size for z counting\n");
  exit(-1);
}

int main (int argc, char *argv[])
{
  double H = 0.50;
  unsigned long int rnd_seed = 0;
  double BASEPOW = DEFAULT_BASEPOW;
  int LAYERSTART = DEFAULT_LAYERSTART;
  int LAYEREND = DEFAULT_LAYEREND;
  const gsl_rng_type * rng_type = gsl_rng_default;
  const char *outfn = DEFAULT_OUTPUT;
  double fieldsize = (1<<20);
  double interval = 1.0;
  unsigned long samplecount = 16;
  double counter_step = 0.50;

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' || argv[i][1] == 0 || argv[i][2] != 0 )
          usage (argv[0]);

      switch (argv[i][1])
        {
          case 'r': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lu", &rnd_seed) != 1 )
                usage (argv[0]);
            break;
          }
          case 'f': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &fieldsize) != 1 )
                usage (argv[0]);
            break;
          }
          case 'i': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &interval) != 1 )
                usage (argv[0]);
            break;
          }
          case 'H': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &H) != 1 )
                usage (argv[0]);
            break;
          }
          case 'c': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &counter_step) != 1 )
                usage (argv[0]);
            break;
          }
          case 's': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lu", &samplecount) != 1 )
                usage (argv[0]);
            break;
          }
          case 'l': {
            if ( i++ == argc || sscanf (argv[i], "%lf", &BASEPOW) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYERSTART) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYEREND) != 1 )
                usage (argv[0]);
            break;
          }
          case 'g': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rng_type = 0;
            for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
              {
                if ( strcmp ((*t)->name, argv[i]) == 0 )
                  {
                    rng_type = *t;
                    break;
                  }
              }
            if ( !rng_type )
              {
                printf ("Random generators available:\n");
                for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
                    printf ("  %s\n", (*t)->name);
                exit(0);
              }
            break;
          }
          case 'o': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            outfn = argv[i];
            break;
          }
          default: {
            usage (argv[0]);
            break;
          }
        }
    }

  assert (LAYERSTART <= LAYEREND);

  rng.init (rng_type, rnd_seed);
  IncrementField field(counter_step, H, BASEPOW, LAYERSTART, LAYEREND);
  zcounter zc(field, fieldsize / counter_step);
  zsampler sampler(zc, interval);
  for ( size_t i = 0; i < samplecount; i++ )
    {
      printf ("zc.sample\n");
      zc.sample();
      printf ("sampler.sample\n");
      sampler.sample(zc.get_range().first, zc.get_range().second);
      printf ("sampler.print\n");
      sampler.print(outfn);
      if ( i != samplecount - 1 )
        {
          printf ("field.reset\n");
          field.reset();
        }
    }

  return 0;
}

